{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Scientific Computing in Python with Numpy\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pre/Post-test\n",
"This test is for testing your current skills in Python. You can use it in two ways:\n",
"* pre-test: to test your skills beforehand. If you are already proficient in Python, and can do this test within approximately 15 minutes, you can scan through the notebook rather than carefully reading each sentence.\n",
"* post-test: to test your skills after Notebook 4. Check whether you learned enough.\n",
"\n",
"**Calculating a derivative**\n",
"Eric is asked to develop a tool for people that need to do some mathematical calculations. He is asked to write a function that calculates and plots the derivative of a given function. The derivative of a continuous function $f(x)$ can be approximated by $f'(x) = \\frac{f(x+\\epsilon)-f(x-\\epsilon)}{2\\epsilon}$ for some small value of $\\epsilon$. As Eric knows that the derivative of $f(x) = \\sin(x)$ is $f'(x) = \\cos(x)$, he uses this function to test whether his function works correct.\n",
"- Make an array in the domain [0, 2*$\\pi$] with 1e4 even spaced values.\n",
"- Plot the function $f(x)$ for this domain.\n",
"- Write a function that calculates the derivative using the approach above.\n",
"- Plot the graph of the derivative in the same figure.\n",
"- Test the correctness of the function by using any other input function.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"### Your Code"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-33372ef05bccfcb0",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Learning objectives\n",
"\n",
"Numpy (Numerical Python) is a library designed for performing scientific computing in Python. \n",
"\n",
"In this notebook, we will introduce numpy arrays, a data structure introduced in numpy for working with vectors and matrices. We will explore how to create them, how to manipulate them, and how to use them for efficient numerical calculations using numpy functions. \n",
"\n",
"After completing this notebook, you are able to:\n",
"* create (multidimensional) numpy arrays from a list of numbers\n",
"* use indexing and slicing with (multidimensional) numpy arrays\n",
"* iterate over a numpy array \n",
"* perform mathematical operations on numpy arrays\n",
"* use functions for creating arrays (eg. `np.zeros()`, `np.linspace()`, `np.random.random()`\n",
"* use numpy functions for vectorized calculations\n",
"* to demonstrate the speed increase of vectorized calculations using `time()`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-7d7b639cd03ba916",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numpy Arrays\n",
"\n",
"You have encountered numpy arrays in previous notebooks. If you don't remember, open the first notebook and read the section where we introduce numpy arrays.\n",
"\n",
"To use numpy arrays, we first need to import the numpy library, which we will do using the shortened name \"np\":"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-4cd897115e9a55da",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2b8a256ccf823e14",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Now that we have imported numpy, we can use functions in numpy to create a numpy array. A simple way to do this is to use the function `np.array()` to make a numpy array from a comma-separated list of numbers in square brackets:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5])\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-cfdb9edf4b2e1707",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Note that numpy does not make a distinction between row vectors and column vectors: they are just vectors. \n",
"\n",
"Look at the cell below, what is the difference with the cell above?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5],dtype=float)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-aa801acee39c518f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Indexing arrays (and counting from zero)\n",
"One useful thing about arrays is that you can access the elements of the array using square brackets:\n",
"\n",
"`a[n]` will give you the n-th element of the array `a`. \n",
"\n",
"This process of extracting a single element from the array is called **indexing**. \n",
"\n",
"Note that here we encounter for the first time what is known as the python **counting from zero** convention. What is the counting from zero convention? In the example above, we created an array:\n",
"\n",
"```\n",
"a = np.array([1,2,3,4,5])\n",
"```\n",
"\n",
"The first element of `a` is `1`. You might think that if you want to access the first element of `a`, you would use the notation `a[1]`. Right?\n",
"\n",
"Let's try it:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[1])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-992a7aa5df6186e4",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**WRONG!** Why? Because the makers of Python decided to start counting from zero: the first element of a tuple `a` is actually `a[0]`. \n",
"\n",
"(This is a long-standing discussion among computer scientists, and the convention is different in many different languages. There are advantages and disadvantages of both, and even essays written about it...but in any case, Python chose to start arrays at zero.)\n",
"\n",
"This also helps better understand the `range()` function: for example, to loop over all the elements in `a`, I can use this code:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(len(a)):\n",
" n = a[i]\n",
" print('a[%d] is %d' % (i,n))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ab7bf8fa5a2ac10f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Here the `len` function returns the length of the array `a`. As we saw before, Python has very smart `for` loops that can automatically iterate over many types of objects, which means we can also print out all the elements of our array like this:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for n in a:\n",
" print(n)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-e1299e2159378fdb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"In Python, if you try to index beyond the end of the array, you will get an error: "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a[5]"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-77aff5e6ccdff810",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"(Remember: indexing starts at zero!)\n",
"\n",
"Python also has a handy feature: negative indices count backwards from the end, and index `-1` corresponds to the last element in the array! "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a[-1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a[-2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8398fd216f6b2f58",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We can also use indexing to change the values of elements in our array:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a)\n",
"a[2] = -1\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3dae75f8494e44f3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.1`** \\\n",
"Set the first three, and the last two, entries of the following array to zero:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5,6,7,8,9,10,11,32,55,78,22,99,55,33.2,55.77,99,101.3])\n",
"\n",
"#some code to set the first three and last two entries to zero \n",
"\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3b080cba952a146f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Slicing numpy arrays\n",
"\n",
"Numpy arrays also support a special type of indexing called \"slicing\" that does not just return a single element of an array, but instead returns a whole part of array. \n",
"\n",
"To do this, I put not just a single number inside my square brackets, but instead two numbers, separated by a colon `:`\n",
"\n",
"`a[n:m]` will return a new tuple that consist of all the elements in `a`, starting at element `n` and ending at element `m-1`. \n",
"\n",
"Let's look at a concrete example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array(range(10))\n",
"print(a)\n",
"print(a[0:5])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-70abffff9b88cf27",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"The notation `a[0:5]` has \"sliced\" out the first five elements of the array.\n",
"\n",
"With slicing, you can also leave off either `n` or `m` from the slice: if leave off `n` it will default to `n=0`, and if you leave off `m`, it will default to the end of the array (also the same as `m=-1` in Python indexing):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[:5])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[5:])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a8b2ffcd10cfc7b1",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Also handy: you can can have Python slice an array with a \"step\" size that is more than one by adding another `:` and a number after that. Find out its operation using:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[0:10:2])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-24b9d17391ceeb7e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Fun: you can also use negative steps:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[-1:-11:-1])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-6ee08109e60cb51f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"And finally, unlike indexing, Python is a bit lenient (merciful) if you slice off the end of an array:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a[0:20])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8d1518c5eafa34e8",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.2`** \\\n",
"Slicing can also be used to *set* multiple values in an array at the same time. Use slicing to set first 10 entries of the array below to zero in one line of code."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array(range(20))+1\n",
"print(a)\n",
"#some code that sets the first 10 entries to zero\n",
"\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**`Exercise 4.2*`** \\\n",
"The same exercise as 4.1 but now code in a smarter way."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#your code here \n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-9b377accb5c0210c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Mathematical operations on arrays\n",
"\n",
"An advantage of using numpy arrays for scientific computing is the way they behave under mathematical operations. In particular, they very often do exactly what we would want them to do if they were a vector:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5])\n",
"print(2*a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a+a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a+1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a-a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a/2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a**2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2963828a5132a24a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"What about if I multiply two vectors together? \n",
"\n",
"In mathematics, if I multiply two vectors, what I get depends on if I use the \"dot product\" or the \"outer product\" for my multiplication:\n",
"\n",
"https://en.wikipedia.org/wiki/Row_and_column_vectors#Operations\n",
"\n",
"The \"dot product\" corresponds to multiplying a column vector by a row vector to produce a single number. The \"outer product\" (also called the \"tensor product\") corresponds to multiplying the column vector by the row vector to make a matrix. \n",
"\n",
"**Question:** If I type `a*a`, or more generally `a*b`, does Python use the inner or outer product? \n",
"\n",
"It turns out: it uses **neither!** In Python, the notation `a*a` produces what is commonly called the \"element-wise\" product: specifically,\n",
"\n",
"`a*b = [a[0]*b[0] a[1]*b[1] a[2]*b[2] ...]`\n",
"\n",
"(Mathematically, this has a fancy name called the Hadamard product, but as you can see, despite the fancy name, it's actually very simple...)\n",
"\n",
"We can see this in action here:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a*a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-713408533e01c84c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"What if I actually want the dot product or the outer product? For that, Python has functions `np.dot()` and `np.outer()`: "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(np.dot(a,a))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(np.outer(a,a))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A useful mathematical operator is the cross product (https://en.wikipedia.org/wiki/Cross_product) where one calculates a vector which is perpendicular to vectors **x** and **y**. \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = np.array([1,0,0])\n",
"y = np.array([0,1,0])\n",
"z = np.cross(x,y)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d29f4ba63452bfee",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Pretty much all operators work with numpy arrays, even comparison operators, which can sometimes be very handy:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(a)\n",
"print(a>3)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8657218bbd2eabbc",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.3`** \\\n",
"Generate a sequence of the first 20 powers of 2 in a numpy array (starting at $2^0$). \n",
"\n",
"Your output should be an array $[2^0, 2^1, 2^2, 2^3, ...]$. \n",
"\n",
"*(Hint: Start with a numpy array created using an appropriate range function that makes an array [0,1,2,3,...])*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# your code that makes the desired array\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-809eedfdadcd8b57",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Functions for creating numpy arrays\n",
"\n",
"In numpy, there are also several handy functions for automatically creating arrays, see https://numpy.org/devdocs/reference/routines.array-creation.html. We provide some examples:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.zeros(5)\n",
"print('zeros', a)\n",
"a = np.ones(5)\n",
"print('ones', a)\n",
"a = np.eye(5)\n",
"print('diagonal')\n",
"print(a)\n",
"a = np.tri(5)\n",
"print('ones below and on diagonal')\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a3a728bdca7ae5e3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### np.linspace\n",
"\n",
"To automatically generate an array with linearly increasing values you can use `np.linspace()`:\n",
"\n",
"https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html\n",
"\n",
"np.linspace takes three arguments: the starting number, the ending number, and the number of points.\n",
"\n",
"This is a bit like the `range` function we saw before, but allows you to pick the total number of points, automatically calculating the (non-integer) step size you need:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.linspace(0,20,40)\n",
"print(a)\n",
"print()\n",
"print(\"Length is: \", len(a))\n",
"print(\"Step size is: \", a[1]-a[0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ea30686a4490a44d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Note that if we wanted to have a step size of exactly 0.5, we need a total of 41 points:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.linspace(0,20,41)\n",
"print(a)\n",
"print()\n",
"print(\"Length is: \", len(a))\n",
"print(\"Step size is: \", a[1]-a[0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-60d973103e7a78a9",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.4`** \\\n",
"Generate an array that runs from -2 to 1 with 20 points using `linspace`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = your code\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-7c4a427b963a9b3c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### np.arange()\n",
"\n",
"If we want to have more control on the exact spacing, we can use the `np.arange()` function. It is like `range()`, asking you for the start, stop, and step size:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.arange(0,20,0.5)\n",
"print(a)\n",
"print()\n",
"print(\"Length is: \", len(a))\n",
"print(\"Step size is: \", a[1]-a[0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5b97dd88408ae9d0",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Here, we already see a small quirk of `arange`: the stopcondition. It does not include this specified number (rather < than <=). If we want to get a range that stops at `20.0`, we need to make the stop point any number a bit bigger than 20 (but smaller than our step size):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.arange(0,20.00000001,0.5)\n",
"print(a)\n",
"print()\n",
"print(\"Length is: \", len(a))\n",
"print(\"Step size is: \", a[1]-a[0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ab1a61494f1f5b33",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"For this reason, I do not find myself using `np.arange()` very often, and mostly use `np.linspace()`. There are also several other useful functions, such as np.geomspace(), which produces geometrically spaced points (such that they are evenly spaced on a log scale). "
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3a0d9465353f954b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.5`** \\\n",
"Generate a numpy array that has a first element with value 60 and last element 50 and takes steps of -0.5 between the values. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = your code\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ed2f97bb11cf695d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Random numbers\n",
"\n",
"Numpy can also generate arrays of random numbers:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.random.random(40)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-624c748e1b750b24",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"This will generate uniform random numbers on the range of 0 to 1, but there are also several other random number generator functions that can make normally distributed random numbers, or random integers, and more."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-bf4eeb182ddf3fdf",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.6`** \\\n",
"Generate a numpy array, `rounded_grades` that contains 300 random grades that have a distribution of a bell-shaped curve that might represent the final grades of the students in this course, with an average grade of 7.5 and a standard deviation of 1. Make sure your grades are rounded to a half point.\n",
"\n",
"(You may find it useful have to look at the help of the `np.random.normal()` function.\n",
"\n",
"(Because of the properties of a normal distribution a small pecentage of the grades may be above a 10, you may leave this for now.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#Your code here that results in a numpy array rounded_grades\n",
"...some code...\n",
"print(rounded_grades)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**`Exercise 4.7`** \\\n",
"There are various ways in which you can analyse your grade distribution.\n",
"\n",
"You can plot a histogram of the grade distribution. Make sure that the width of our histogram bars is 0.5, corresponding to our rounded grades."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"... some code ..."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-042c8d00bb09795d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Multidimensional arrays (matrices)\n",
"\n",
"\n",
"We have looked at 1D arrays especially. However, numpy also supports two-dimensional (or N-dimensional) numpy arrays, that can represent matrices. To make a 2D numpy array, you can use the `zeros()` function, for example, but with a two-entry list of numbers specifying the size N and M of the matrix:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = np.zeros([10,10])\n",
"print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-63e005300e7f7098",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"For two dimensional matrices, the usual function `len()` is not enough to tell us about the shape of our matrix. Instead, we can use a property of the numpy matrix itself called its `shape`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(len(m))\n",
"print(m.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f8e1398b61d09699",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Indexing two dimensional arrays works by using commas inside square brackets to specify the index of the first and second dimensions:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array(range(1,6))\n",
"m = np.outer(a,a)\n",
"print(\"Initial matrix:\")\n",
"print(m)\n",
"\n",
"# First index in the row number (counting from zero), second index in the column number\n",
"m[2,3] = -1 \n",
"print(\"\\nAfter changing entry [2,3] to -1:\")\n",
"print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-4e9af994bf9d7df7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"You can also use slicing to to assign values to an array from a vector, which can be a handy way to enter a matrix by hand:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = np.zeros([3,3])\n",
"m[0,:] = [1,2,3]\n",
"m[1,:] = [4,5,6]\n",
"m[2,:] = [7,8,9]\n",
"print(m)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-461ca601abc17699",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Similarly, slicing also can be used to extract rows, columns, or blocks of the matrix:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# A row\n",
"print(m[1,:])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# A column\n",
"print(m[:,1])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Extract a block as a sub-matrix\n",
"print(m[1:,1:])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-6ab7d9e10575586d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"There are several functions for making matrices which you may find useful someday: \n",
"\n",
"https://docs.scipy.org/doc/numpy/reference/routines.array-creation.html\n",
"\n",
"including this one which is used often:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# The identity matrix\n",
"print(np.eye(10))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# A band diagonal matrix\n",
"print(np.eye(10,k=-1))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-484f34f51c8f8c5e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.8`** \\\n",
"Use Python to calculate the following matrix multiplication: \n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
"1 & 1 & 0 \\\\\n",
"0 & 2 & 1 \\\\\n",
"1 & 0 & 1\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
"1 & 3 & 0 \\\\\n",
"3 & 1 & 1 \\\\\n",
"1 & 1 & 1\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"To perform this multiplication, you will need to use the `matmul` or `dot` routine of numpy: \n",
"\n",
"https://docs.scipy.org/doc/numpy/reference/generated/numpy.matmul.html\n",
"https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html\n",
"\n",
"*(For the nerds: do these two matrices commute?)*\n",
"\n",
"*(For the real nerds: have your program check if they commute and inform the user!)*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m1 = np.zeros([3,3])\n",
"some code to fill the matrix\n",
"m2 = np.zeros([3,3])\n",
"some code to fill the matrix\n",
"\n",
"# Check the matrices \n",
"print(m1)\n",
"print()\n",
"print(m2)\n",
"print()\n",
"\n",
"product = some code\n",
"print(product)\n",
"\n",
"do nerd stuff if you want..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can do seemingly smart things using these 'special' vectors and matrices. For instance, what if we want to take the sum of all values in an array? We can use the dotproduct:\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
"1 & 1 & 1 \\\\\n",
"\\end{bmatrix}\n",
"\\cdot\n",
"\\begin{bmatrix}\n",
"1 \\\\ 2 \\\\ 3 \\\\\n",
"\\end{bmatrix}\n",
"= 6\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ones = np.ones(3)\n",
"array = np.linspace(1,3,3)\n",
"print(np.dot(one_diag,array))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover, if we want to create a cummulative sum of values in an array, we can use the tri matrix:\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
"1 & 0 & 0 \\\\\n",
"1 & 1 & 0 \\\\\n",
"1 & 1 & 1\n",
"\\end{bmatrix}\n",
"\\cdot\n",
"\\begin{bmatrix}\n",
"1 \\\\\n",
"2 \\\\\n",
"3 \n",
"\\end{bmatrix}\n",
"= \n",
"\\begin{bmatrix}\n",
"1\\cdot1 + 0\\cdot2 + 0\\cdot3 \\\\\n",
"1\\cdot1 + 1\\cdot2 + 0\\cdot3 \\\\\n",
"1\\cdot1 + 1\\cdot2 + 1\\cdot3\n",
"\\end{bmatrix}\n",
"=\n",
"\\begin{bmatrix}\n",
"1 \\\\\n",
"3 \\\\\n",
"6 \n",
"\\end{bmatrix}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tri = np.tri(3)\n",
"print(np.dot(tri,array))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another way to do so is to use a for-loop: "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.linspace(1,3,3)\n",
"sum_a = 0\n",
"cumsum = np.array([])\n",
"\n",
"for x in a:\n",
" sum_a += x\n",
" cumsum = np.append(cumsum,sum_a)\n",
" \n",
"print(\"The sum is\", sum_a)\n",
"print(\"The cummulative sum is\", cumsum)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**`Exercise 4.9`** \n",
"\n",
"Below we have a set of repeated measurements. We want to understand how the average value of our repeated measurements evaluates over time, that is, how does the average value change as function of $N$ repeated measurements? To investigate this: calculate and plot the average value as function of the number of measurements."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"measurements = np.random.normal(5,1,20)\n",
"print(measurements)\n",
"\n",
"#your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But wait... is there not a more direct way to do all of the above... Sure! We can make use of Numpy functions."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0a8f096130330d1f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numpy functions\n",
"\n",
"For many common (mathematical) operations, fuctions exists. For instance, we can calculate the average value of our measurements of above using a for loop, but also using the function `np.average`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sum_meas = 0\n",
"\n",
"for measurement in measurements:\n",
" sum_meas += measurement\n",
"\n",
"print(sum_meas/len(measurements))\n",
"print(np.average(measurements))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-48cd73f9aeaf481e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"This is very handy: it saves us loads of thinking and typing! From the function name, it is also easy to understand what you are doing, making the code clearer and easier to read. However, the purpose of numpy functions is not only to save lots of typing: they also can often perform calculations MUCH faster than if you do program the calculation yourself with a `for` loop, as we will see in the next section.\n",
"\n",
"Python also has many other useful functions for performing calculations using arrays:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Calculate the standard deviation\n",
"print(np.std(measurements,ddof=1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# The square root\n",
"print(np.sqrt(a))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Numpy also has max and min, here is an example of min\n",
"a = np.linspace(-10,10,100)\n",
"print(np.min(a**2))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3c32880c009136b2",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Good question for you to think about: why is the minimum value not zero? And what would I have to change above to get the code to return zero? \n",
"\n",
"In addition to finding the minimum value in a vector, the function `argmin` can tell you **where** (what index number) the minimum is:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Find the index number of the minimum of the array\n",
"i = np.argmin(a**2)\n",
"print(i)\n",
"print((a**2)[i])"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d520c5750ce4f54a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Note also here that we used round brackets `()` around the `a**2` in the `print` statement to be able to then index the resulting array `a**2` array using square brackets `[]`. \n",
"\n",
"You can find the full list of mathematical numpy functions on the documentation website:\n",
"\n",
"https://docs.scipy.org/doc/numpy/reference/routines.math.html\n",
"\n",
"and the full list of all functions in the reference guide: \n",
"\n",
"https://docs.scipy.org/doc/numpy/reference/index.html"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-081110b083e1eeb8",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**`Exercise 4.10`** \n",
"\n",
"**a)**\n",
"Create an array with values ranging from 14 to 42, skipping the odd numbers.\n",
"\n",
"**b)**\n",
"Make an array `x` that runs from 0 to 4 with 20 points, and calculate an array `y` whose entries are equal to the square root of the entries in `x`. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"your code to make the requested arrays x and y\n",
"print(y)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-dfed25e05b150ef3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"### Vectorisation\" and fast code with numpy functions\n",
"\n",
"\n",
"In the first example above, we showed two ways of calculating an average: one using a `for` loop, and one using the numpy function. \n",
"\n",
"Functionally, they are equivalent: they do exactly the same thing. \n",
"\n",
"A curious feature of Python is that if you use functions instead of coding loops yourself, often things are **MUCH MUCH** faster. \n",
"\n",
"To show this quantitatively, we will use the `time` library to calculate the time it takes to find the average of a pretty big array using both techniques:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# The time() function from the time library will return a floating point number representing \n",
"# the number of seconds since January 1, 1970, 00:00:00 (UTC), with millisecond or even microsecond\n",
"# precision\n",
"# \n",
"# We will use this to make a note of the starting time and the ending time, \n",
"# and then print out the time difference \n",
"from time import time\n",
"\n",
"# A pretty big array, 50 million random numbers\n",
"a = np.random.random(int(50e6))\n",
"\n",
"# Set timer\n",
"t1 = time()\n",
"\n",
"# Caluclate the average\n",
"avg = 0\n",
"for x in a:\n",
" avg += x \n",
"avg /= len(a)\n",
"t2 = time()\n",
"t_forloop = t2-t1\n",
"print(\"The 'for' loop took %.3f seconds\" % (t2-t1))\n",
"\n",
"t1 = time()\n",
"avg = np.average(a)\n",
"t2 = time()\n",
"t_np = t2-t1\n",
"print(\"Numpy took %.3f seconds\" % (t_np))\n",
"\n",
"# Now let's compare them\n",
"print(\"\\nNumpy was %.1f times faster!\" % (t_forloop/t_np))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-bbb6d85899ba86b1",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Why is numpy so much faster? The reason is that Python is an interpreted language. In each of the steps of the `for` loop, the Python kernel reads in the next step it has to do, translates that into an instruction for your computer processor, asks the computer to perform the step, gets the result back, reads in the next step, translates that into a processor instruction, sends that as an instruction to the computer processor, etc, etc. \n",
"\n",
"If we did the same test in a compiled programing language like C, there would be no difference if we used a library function or if we wrote our own `for` loop. \n",
"\n",
"When you use smart functions in Python libraries, like (many of) those in numpy, numpy will actually use an external library compiled in a language like C or Fortran that is able to send all of the calculation in one step to your computer processor, and in one step, get all the data back. This makes Python nearly as fast as a compiled language like C or Fortran, as long as you are smart in how you use it and avoid having \"manual\" `for` loops for large or long calculations. \n",
"\n",
"(For small calculations, Python coded `for` loops are perfectly fine and very handy!)\n",
"\n",
"In the language of interpreted programmers, finding smart ways of getting what you need done using \"compiled library functions\" is often referred to as vectorisation. \n",
"\n",
"Note that even normal mathematical operators are actually \"vectorized functions\" when they operate:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# This is actually a vectorized 'for' loop, it involves multiplying 50 million numbers by 5\n",
"b = 5*a"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2a60c21baccb3c7d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
" Here is a nice example of a vectorized way of counting the number of times the number '5' occurs in a random sample of 100 integers between 0 and 20:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nums = np.random.randint(0,21,100)\n",
"print(\"There are %d fives\" % np.sum(nums == 5))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c8bc80849f5e0e0a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"To see how this works, we can look at the intermediate steps:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nums = np.random.randint(0,21,100)\n",
"print(nums)\n",
"print(nums == 5)\n",
"print(\"There are %d fives\" % np.sum(nums == 5))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-1a3f30004c3e3ffb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Note that in this case, `np.sum()` will convert the `bool` value `True` into `1` and `False` into `0` for calculating the sum, according the the standard convertion of `bool` types to `int` types. You can see this in action if you want using the function `astype()` that is built into numpy arrays:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(nums == 5)\n",
"print((nums == 5).astype('int'))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-dd8dc8bd0ab88a1c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Another neat feature that `numpy` has is that is can 'vectorize' normal Python functions so that they can take `numpy` functions and make them a bit faster (10-20%). This is done using the `np.frompyfunc` function. An example is given below."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def fac(N):\n",
" if (N == 1) or (N == 0):\n",
" return 1\n",
" return N * fac(N-1)\n",
"\n",
"#As you can see applying the following array does not work\n",
"N = np.arange(9)\n",
"\n",
"try:\n",
" print(fac(N))\n",
" print(\"This works\")\n",
"except:\n",
" print(\"This does not work \\n\")\n",
" \n",
"#Now if we make this a numpy function is will work \n",
"#(the numbers in the frompyfunc represent the input and output values)\n",
"\n",
"numpy_fac = np.frompyfunc(fac,1,1)\n",
"\n",
"\n",
"try:\n",
" print(numpy_fac(N))\n",
" print(\"This works\")\n",
"except:\n",
" print(\"This does not work\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Monte Carlo Simulation\n",
"\n",
"\n",
"A commonly used method for numerical calculations is a Monte Carlo simulation (named after the (in)famous casino district). Simply stated, the method consists of generating a great number of random points and checking afterwards which points satisfy the boundary conditions of the calculation.\n",
"\n",
"For example: to calculate the area of a circle with radius one, one can generate a great number of random points within a square of size 1 by 1, and check for each points whether they fall within the circle $r < 1$. The area of the circle can then be obtained by multiplying the area of the square with the number of points that are within the circle divided by the total number of points.\n",
"\n",
"This may sound a bit cumbersome, and of course in this example calculating $\\pi \\cdot r^2$ is a lot faster. There are however many applications where the Monte Carlo method proves to be very useful; for example in calculating integrals that cannot (easily) be solved analytically.\n",
"\n",
"In the next exercises, you are going to calculate the integral of the function $f(x) = e^{x^2}$ in the interval [0,1] using the Monte Carlo method. \n",
"\n",
"**Exercise 4.11** \n",
"\n",
"First, make a plot to get an idea of what the function looks like. \n",
"\n",
"What is the maximum value of $f(x)$ in the interval [0,1]?\n",
"\n",
"Make a (rough) estimate of the value of the integral by looking a the graph"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"#First we define our function\n",
"def f(x):\n",
" return np.exp(x**2)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are going to generate random samples in a square that has an area that is greater than the integral that we want to calculate. \n",
"\n",
"For this, we will use the np.random.uniform() function, which generates random float in the interval [a,b) - notice, this is an half-open interval so b is not included. \n",
"\n",
"Moreover, it is wise to define the area by four points: $x_a$, $x_b$, $y_a$ and $y_b$.\n",
"\n",
"**Exercise 4.12**\n",
"\n",
"Complete the code below to generate a random sample of y-values in an appropriate range."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"N = int(1e6) #number of points that we are going to use in our calculation\n",
"x_a = \n",
"x_b = \n",
"y_a = \n",
"y_b = \n",
"#Generate samples:\n",
"random_x = np.random.uniform( , ,size=N) \n",
"random_y = np.random.uniform( , ,size=N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we want to determine the number of points in our sample that fall within the area under the line f(x). \n",
"\n",
"**Exercise 4.13**\n",
"\n",
"Determine the number of points that satisfy this condition using a for loop and calculate the value of the integral."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.464291\n",
"The solution of the integral is 1.464291\n",
"Time for calculation: 5.187 s\n"
]
}
],
"source": [
"from time import time\n",
"\n",
"t1 = time()\n",
"\n",
"### your code ###\n",
"\n",
"\n",
"\n",
"solution_integral = \n",
"\n",
"print('The solution of the integral is %.6f' %(solution_integral)) \n",
"print('Time for calculation: %.3f s' %(time()-t1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have a solution of the integral (much easier than on paper, right?)! But how do we know how accurate our solution is? First, compare your result to your estimate from the previous exercise. Does it make sense?\n",
"\n",
"To make a substantiated claim about the accuracy of our solution we have to determine the uncertainty in our result. As we are basically performing a count, Poisson statistics tells us that the uncertainty in the counted number of points that satisfy the condition is the square root of the counted number: $u(N_{counted}) = \\sqrt{N_{counted}}$.\n",
"We can then determine the uncertainty in the area using the calculus approach.\n",
"\n",
"Let's check our solution by comparing it to the answer given by another numerical method from the scipy library."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Our solution: 1.464 +/- 0.002\n",
"Scipy solution: 1.46265174590718 +/- 0.00000000000002\n",
"Difference between the solutions: 0.001639\n"
]
}
],
"source": [
"from scipy.integrate import quad\n",
"\n",
"Area = np.abs(x_b-x_a)*np.abs(y_b-y_a)\n",
"err_solution = Area*np.sqrt(s)/N #estimated error of our solution, which is the area of our random uniform box divided by the number of points used in the calculation times the poisson uncertainty\n",
"\n",
"scipy_solution = quad(f,0,1) #calculating integral using quad function from scipy, returns: value of integral and estimated error\n",
"\n",
"abs_difference = np.abs(solution_integral-scipy_solution[0]) #determine absolute difference between solutions\n",
"\n",
"\n",
"assert abs_difference <= 2*np.sqrt(err_solution**2 + scipy_solution[1]**2), 'The results are not in agreement' #check if the values are in agreement \n",
"\n",
"print('Our solution: %.3f +/- %.3f' %(solution_integral, err_solution))\n",
"print('Scipy solution: %.14f +/- %.14f' %(scipy_solution[0], scipy_solution[1]))\n",
"print('Difference between the solutions: %.6f' %(abs_difference))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looks pretty good!\n",
"\n",
"Our method gives a reasonable answer, but as you can see the uncertainty is relatively high...If we want to improve our result we should use more points to lower the uncertainty (which scales with $\\sqrt{N}$). However, the for loop we used is already quite slow. More points would mean an annoyingly long calculation. \n",
"\n",
"We should therefore try to avoid the use of a for loop and instead use a comparison operator to speed up the calculation.\n",
"\n",
"**Exercise 4.14**\n",
"\n",
"Again, complete the code below and calculate the integral, this time without using a for loop"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t1 = time()\n",
"\n",
"s = np.sum( <= )\n",
"\n",
"solution_integral = \n",
"\n",
"print('The solution of the integral is %.6f' %(solution_integral)) \n",
"print('Time for calculation: %.3f s' %(time()-t1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is much faster!\n",
"\n",
"Now lastly we want to plot our random generated points that are within the area underneath f(x) and those that are not in the same graph with a different colour. For this, we need the actual values of the points that satisfy the condition.\n",
"\n",
"As you have seen before, a comparison operation on a numpy array returns a type boolean array:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"array = np.array([4,2,6,8,5,4])\n",
"\n",
"boolean_array = array > 4\n",
"\n",
"print(boolean_array)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can subsequently use this boolean array to get the values of the array that satisfy the condition:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"new_array = array[boolean_array]\n",
"\n",
"print(new_array)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or, more simply:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"new_array = array[array > 4]\n",
"\n",
"print(new_array)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.15**\n",
"\n",
"Make a figure in which you plot:\n",
"\n",
"- the function f(x)\n",
"- the random points that are within the area under the line f(x)\n",
"- the random points that are outside this area (in a different colour)\n",
"\n",
"We will generate new random samples with fewer points, otherwise we cannot distinguish them.\n",
"\n",
"Don't forget to label your axes!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"random_x = np.random.uniform( , ,1000)\n",
"random_y = np.random.uniform( , ,1000)\n",
"\n",
"#your code"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.16** \n",
"\n",
"One way to 'improve' the quality of your data, making it less noisy, is the use of a moving average filter. The moving average filter calculates the average value of N elements, N often being 3 or 5. This smooth out any abrupt changes and allows a better focus on long term trends. Of course you will lose two elements of the entire dataset (think yourselves why). \n",
"\n",
"\n",
"The filtered data is stored in an array: $F(i) = \\frac{p(i)+p(i+1)+p(i+2)}{3}$, looped over all elements $p(i)$. More precise: $F(i)=\\frac{1}{k}\\sum_{i}^{i+k-1}p(i)$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Z = np.arange(1,20)\n",
"#def moving_average_1(p, k):\n",
"# F = np.empty(len(p)-(k-1))\n",
"\n",
"print(Z)\n",
"print(moving_average_1(Z,3))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2f8479d8cc410928",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Solutions to exercises"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0a4e3eec30674d19",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.1** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-daf14d85fa463aba",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5,6,7,8,9,10,11,32,55,78,22,99,55,33.2,55.77,99,101.3])\n",
"\n",
"#some code to set the first three and last two entries to zero \n",
"\n",
"# A bit tricky, but these work\n",
"a[0] = 0\n",
"a[1] = 0\n",
"a[2] = 0\n",
"\n",
"# We could count forwards (a[19] = 0 and a[20] = 0), but is much \n",
"# easier to count backwards from the end\n",
"a[-2] = 0\n",
"a[-1] = 0\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-1f286c69f57cbc80",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.2:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-efb5e7292fc18003",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a = np.array(range(20))+1\n",
"print(a)\n",
"a[0:10] = 0\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.2$*$:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([1,2,3,4,5,6,7,8,9,10,11,32,55,78,22,99,55,33.2,55.77,99,101.3])\n",
"a[-2:]=0\n",
"a[:2]=0\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b47cf723add61750",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.3:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-87b5343c8e0fdbce",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"n = np.array(range(21))\n",
"out = 2**n\n",
"print(out)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ea7ddbad81cf4ad1",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.4:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-e6d951d0073dc45f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a = np.linspace(-2,1,20)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d33f7386d70b8e4f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.5:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-84bb816a29f6f1bb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"a = np.arange(60,49.9,-0.5)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-095bdb4fd5d1fae4",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.6:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5a51a5dde00116ba",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"raw = np.random.normal(7.5,1,300)\n",
"rounded_grades = np.round(raw*2)/2\n",
"print(rounded_grades)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f8253b40b7bb046e",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.7:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"bins = np.arange(np.min(rounded_grades),np.max(rounded_grades)+0.5,0.5)\n",
"\n",
"plt.figure()\n",
"plt.hist(rounded_grades, bins=bins)\n",
"plt.xlabel('Grade')\n",
"plt.ylabel('#')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.8:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a69f7838d0b9e25a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"m1 = np.zeros([3,3])\n",
"m1[0,:] = np.array([1,1,0])\n",
"m1[1,:] = np.array([0,2,1]) \n",
"m1[2,:] = np.array([1,0,1]) \n",
"m2 = np.zeros([3,3])\n",
"m2[:,0] = np.array([1,3,1])\n",
"m2[:,1] = np.array([3,1,1]) \n",
"m2[:,2] = np.array([0,1,1]) \n",
"\n",
"# Check the matrices \n",
"print(m1)\n",
"print()\n",
"print(m2)\n",
"print()\n",
"\n",
"product = np.matmul(m1,m2)\n",
"print(product)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.9:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"measurements = np.random.normal(5,1,20)\n",
"print(measurements)\n",
"\n",
"average_n = np.cumsum(measurements)/np.arange(1,len(measurements)+1)\n",
"\n",
"plt.figure()\n",
"plt.plot(np.arange(1,len(measurements)+1), average_n)\n",
"plt.xlabel('average')\n",
"plt.ylabel('n')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3ba058efb16a91ad",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"**Exercise 4.10:** "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-57c551e2a58429d8",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"array_skip = np.arange(14,44,2)\n",
"print(array_skip)\n",
"\n",
"\n",
"x = np.linspace(0,4,20)\n",
"y = np.sqrt(x)\n",
"print(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Excercise 4.11**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"#First we define our function\n",
"def f(x):\n",
" return np.exp(x**2)\n",
"\n",
"x = np.linspace(0,1,1000) #generate an array with x values\n",
"\n",
"#make the figure\n",
"plt.figure(figsize=(6,4))\n",
"plt.plot(x,f(x))\n",
"plt.xlim(0,1)\n",
"plt.ylim(0)\n",
"plt.ylabel('f(x)')\n",
"plt.xlabel('x')\n",
"plt.show()\n",
"\n",
"print('Maximum value of f(x) in the interval [0,1] is %.3f' %(np.max(f(x))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The area under the graph is roughly equal to half of the total area of the square 1x3, i.e. 1.5."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.12**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"N = int(1e6) #number of points that we are going to use in our calculation\n",
"x_a = 0\n",
"x_b = 1\n",
"y_a = 0\n",
"y_b = 3\n",
"#Generating our samples:\n",
"random_x = np.random.uniform(x_a,x_b,size=N) \n",
"random_y = np.random.uniform(y_a,y_b,size=N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.13**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from time import time\n",
"\n",
"t1 = time()\n",
"\n",
"s = 0\n",
"\n",
"for i in range(len(random_x)):\n",
" if random_y[i] <= np.exp(random_x[i]**2):\n",
" s+= 1\n",
"\n",
"solution_integral = (x_b-x_a)*(y_b-y_a)*s/N\n",
"\n",
"print('The solution of the integral is %.6f' %(solution_integral)) \n",
"print('Time for calculation: %.3f s' %(time()-t1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.14**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t1 = time()\n",
"\n",
"s = np.sum(random_y <= np.exp(random_x**2))\n",
"\n",
"solution_integral = 3*s/N\n",
"\n",
"print('The solution of the integral is %.6f' %(solution_integral)) \n",
"print('Time for calculation: %.3f s' %(time()-t1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.15**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"random_x = np.random.uniform(0,1,1000)\n",
"random_y = np.random.uniform(0,3,1000)\n",
"\n",
"points_within_area = np.array([random_x[random_y <= np.exp(random_x**2)],random_y[random_y <= np.exp(random_x**2)]])\n",
"points_outside_area = np.array([random_x[random_y > np.exp(random_x**2)],random_y[random_y > np.exp(random_x**2)]])\n",
"\n",
"plt.figure(figsize=(6,4))\n",
"plt.plot(points_outside_area[0],points_outside_area[1], '.',c='red')\n",
"plt.plot(points_within_area[0],points_within_area[1], '.', c='green')\n",
"plt.plot(x,f(x),linewidth=3,label='f(x)')\n",
"plt.xlim(0,1)\n",
"plt.ylim(0,3)\n",
"plt.ylabel('f(x)')\n",
"plt.xlabel('x')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise 4.16**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def moving_avg(p, k):\n",
" \n",
" F = np.empty(len(p)-(k-1))\n",
" \n",
" for i in range(len(p)-(k-1)):\n",
" F[i] = (1/k)*np.sum(p[i:i+k])\n",
" \n",
" return F\n",
"\n",
"Z = np.arange(4,20)\n",
"print(Z)\n",
"print(moving_average(Z, 3))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,md"
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.16"
}
},
"nbformat": 4,
"nbformat_minor": 4
}